// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

/* 
 
 * NOTE: This file is the modified version of [s,d,c,z]column_dfs.c file in SuperLU 
 
 * -- SuperLU routine (version 2.0) --
 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
 * and Lawrence Berkeley National Lab.
 * November 15, 1997
 *
 * Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
 *
 * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
 * EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
 *
 * Permission is hereby granted to use or copy this program for any
 * purpose, provided the above notices are retained on all copies.
 * Permission to modify the code and to distribute modified code is
 * granted, provided the above notices are retained, and a notice that
 * the code was modified is included with the above copyright notice.
 */
#ifndef SPARSELU_COLUMN_DFS_H
#define SPARSELU_COLUMN_DFS_H

template <typename Scalar, typename StorageIndex> class SparseLUImpl;
namespace Eigen {

namespace internal {

    template <typename IndexVector, typename ScalarVector> struct column_dfs_traits : no_assignment_operator
    {
        typedef typename ScalarVector::Scalar Scalar;
        typedef typename IndexVector::Scalar StorageIndex;
        column_dfs_traits(Index jcol, Index& jsuper, typename SparseLUImpl<Scalar, StorageIndex>::GlobalLU_t& glu, SparseLUImpl<Scalar, StorageIndex>& luImpl)
            : m_jcol(jcol), m_jsuper_ref(jsuper), m_glu(glu), m_luImpl(luImpl)
        {
        }
        bool update_segrep(Index /*krep*/, Index /*jj*/) { return true; }
        void mem_expand(IndexVector& lsub, Index& nextl, Index chmark)
        {
            if (nextl >= m_glu.nzlmax)
                m_luImpl.memXpand(lsub, m_glu.nzlmax, nextl, LSUB, m_glu.num_expansions);
            if (chmark != (m_jcol - 1))
                m_jsuper_ref = emptyIdxLU;
        }
        enum
        {
            ExpandMem = true
        };

        Index m_jcol;
        Index& m_jsuper_ref;
        typename SparseLUImpl<Scalar, StorageIndex>::GlobalLU_t& m_glu;
        SparseLUImpl<Scalar, StorageIndex>& m_luImpl;
    };

    /**
 * \brief Performs a symbolic factorization on column jcol and decide the supernode boundary
 * 
 * A supernode representative is the last column of a supernode.
 * The nonzeros in U[*,j] are segments that end at supernodes representatives. 
 * The routine returns a list of the supernodal representatives 
 * in topological order of the dfs that generates them. 
 * The location of the first nonzero in each supernodal segment 
 * (supernodal entry location) is also returned. 
 * 
 * \param m number of rows in the matrix
 * \param jcol Current column 
 * \param perm_r Row permutation
 * \param maxsuper  Maximum number of column allowed in a supernode
 * \param [in,out] nseg Number of segments in current U[*,j] - new segments appended
 * \param lsub_col defines the rhs vector to start the dfs
 * \param [in,out] segrep Segment representatives - new segments appended 
 * \param repfnz  First nonzero location in each row
 * \param xprune 
 * \param marker  marker[i] == jj, if i was visited during dfs of current column jj;
 * \param parent
 * \param xplore working array
 * \param glu global LU data 
 * \return 0 success
 *         > 0 number of bytes allocated when run out of space
 * 
 */
    template <typename Scalar, typename StorageIndex>
    Index SparseLUImpl<Scalar, StorageIndex>::column_dfs(const Index m,
                                                         const Index jcol,
                                                         IndexVector& perm_r,
                                                         Index maxsuper,
                                                         Index& nseg,
                                                         BlockIndexVector lsub_col,
                                                         IndexVector& segrep,
                                                         BlockIndexVector repfnz,
                                                         IndexVector& xprune,
                                                         IndexVector& marker,
                                                         IndexVector& parent,
                                                         IndexVector& xplore,
                                                         GlobalLU_t& glu)
    {
        Index jsuper = glu.supno(jcol);
        Index nextl = glu.xlsub(jcol);
        VectorBlock<IndexVector> marker2(marker, 2 * m, m);

        column_dfs_traits<IndexVector, ScalarVector> traits(jcol, jsuper, glu, *this);

        // For each nonzero in A(*,jcol) do dfs
        for (Index k = 0; ((k < m) ? lsub_col[k] != emptyIdxLU : false); k++)
        {
            Index krow = lsub_col(k);
            lsub_col(k) = emptyIdxLU;
            Index kmark = marker2(krow);

            // krow was visited before, go to the next nonz;
            if (kmark == jcol)
                continue;

            dfs_kernel(StorageIndex(jcol), perm_r, nseg, glu.lsub, segrep, repfnz, xprune, marker2, parent, xplore, glu, nextl, krow, traits);
        }  // for each nonzero ...

        Index fsupc;
        StorageIndex nsuper = glu.supno(jcol);
        StorageIndex jcolp1 = StorageIndex(jcol) + 1;
        Index jcolm1 = jcol - 1;

        // check to see if j belongs in the same supernode as j-1
        if (jcol == 0)
        {  // Do nothing for column 0
            nsuper = glu.supno(0) = 0;
        }
        else
        {
            fsupc = glu.xsup(nsuper);
            StorageIndex jptr = glu.xlsub(jcol);  // Not yet compressed
            StorageIndex jm1ptr = glu.xlsub(jcolm1);

            // Use supernodes of type T2 : see SuperLU paper
            if ((nextl - jptr != jptr - jm1ptr - 1))
                jsuper = emptyIdxLU;

            // Make sure the number of columns in a supernode doesn't
            // exceed threshold
            if ((jcol - fsupc) >= maxsuper)
                jsuper = emptyIdxLU;

            /* If jcol starts a new supernode, reclaim storage space in
     * glu.lsub from previous supernode. Note we only store 
     * the subscript set of the first and last columns of 
     * a supernode. (first for num values, last for pruning)
     */
            if (jsuper == emptyIdxLU)
            {  // starts a new supernode
                if ((fsupc < jcolm1 - 1))
                {  // >= 3 columns in nsuper
                    StorageIndex ito = glu.xlsub(fsupc + 1);
                    glu.xlsub(jcolm1) = ito;
                    StorageIndex istop = ito + jptr - jm1ptr;
                    xprune(jcolm1) = istop;  // initialize xprune(jcol-1)
                    glu.xlsub(jcol) = istop;

                    for (StorageIndex ifrom = jm1ptr; ifrom < nextl; ++ifrom, ++ito) glu.lsub(ito) = glu.lsub(ifrom);
                    nextl = ito;  // = istop + length(jcol)
                }
                nsuper++;
                glu.supno(jcol) = nsuper;
            }  // if a new supernode
        }      // end else:  jcol > 0

        // Tidy up the pointers before exit
        glu.xsup(nsuper + 1) = jcolp1;
        glu.supno(jcolp1) = nsuper;
        xprune(jcol) = StorageIndex(nextl);  // Initialize upper bound for pruning
        glu.xlsub(jcolp1) = StorageIndex(nextl);

        return 0;
    }

}  // end namespace internal

}  // end namespace Eigen

#endif
